Access to care

library(tidyverse)
library(rlang)
library(leaflet)

hospitals <- read_csv("../data/hospital_list.csv")

Parsed with column specification:
cols(
  .default = col_character(),
  zip = col_double(),
  `Phone Number` = col_double(),
  score = col_double(),
  longitude = col_double(),
  latitude = col_double(),
  population = col_double(),
  poverty = col_double(),
  percent = col_double(),
  poverty_percent = col_double()
)
See spec(...) for full column specifications.

hospitals

hospital_locations <- hospitals  %>%
  select(state, longitude, latitude)
write_rds(hospital_locations, "../data/hospital_locations.rds")
hospital_locations

population <- usmap::countypop

hospital_count <- hospitals %>%
  count(state, original_county) 


state_hospitals <- population %>%
  left_join(hospital_count, by = c("abbr" = "state", "county" = "original_county")) %>%
  mutate(hospitals = ifelse(is.na(n), 0, n)) %>%
  select(
    fips, 
    state = abbr, 
    county, 
    population = pop_2015,
    hospitals
    )

state_hospitals

set.seed(100)

hospital_sample <- state_hospitals %>%
  sample_frac(0.3)

hospital_sample

model <- lm(hospitals ~ population, data = hospital_sample)
write_rds(model, "../data/model.rds")
summary(model)

Call:
lm(formula = hospitals ~ population, data = hospital_sample)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.7909  -0.7106   0.0908   0.2498   8.3805 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.027e-01  3.942e-02   17.83   <2e-16 ***
population  7.907e-06  1.182e-07   66.90   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.149 on 941 degrees of freedom
Multiple R-squared:  0.8263,    Adjusted R-squared:  0.8261 
F-statistic:  4475 on 1 and 941 DF,  p-value: < 2.2e-16

predictions <- predict(model, 
  newdata = state_hospitals, 
  interval = "prediction") %>%
  as_tibble() %>%
  mutate_all(round)

predictions

hospital_results <- state_hospitals %>%
  bind_cols(predictions) %>%
  mutate(result = case_when(
    hospitals < lwr ~ -1,
    hospitals > upr ~ 1,
    TRUE ~ 0)) %>%
  mutate(
    county_key = str_remove_all(county, " County"),
    county_key = str_remove_all(county_key, " Parish"),
    county_key = str_remove_all(county_key, " city"),
    county_key = str_to_lower(county_key),
    county_key = str_remove_all(county_key, " "),
    county_key = str_replace_all(county_key, "st. ", "saint"),
      )  

write_rds(hospital_results, "../data/hospitals.rds")

hospital_results 

shapes <- read_rds("../data/shapes.rds") %>%
  mutate(
    county_key = str_to_lower(county),
    county_key = str_remove_all(county_key, " "),
    county_key = str_replace_all(county_key, "st. ", "saint")
  ) 

shapes 

county_hospitals <- hospital_results  %>%
  select(- county) %>%
  left_join(shapes, by = c("state", "county_key")) %>%
  filter(state != "PR", state != "AK") %>%
  filter(!is.na(county))

write_rds(county_hospitals, "../data/county_hospitals.rds", compress = "gz")

county_hospitals

states <- county_hospitals %>%
  group_by(state) %>%
  summarise() %>%
  pull()

under <- "#CC79A7"
over <- "#0072B2"
at_level <- "#008b00"
hospital_color <- "#F0E442"
st <- c("LA")
#st <- states
include <- c(-1, 1, 0)
add_hospitals <- 0

initial_map <- leaflet() %>%
  addProviderTiles(providers$CartoDB) %>%
  addLegend("bottomright", 
            color  = c(under,over ,at_level, hospital_color), 
            labels = c("Less hopitals than expected",
                       "More hospitals than expected", "Within Range", 
                       "Hospital Location"), 
            title  = "Legend",opacity = 0.5)

locations <- hospital_locations %>%
  filter(state %in% st) %>%
  select(longitude, latitude) %>%
  mutate_all(round, 3) %>%
  count(longitude, latitude) %>%
  mutate(popup = paste0("Hospitals: ", n))

counties <- county_hospitals %>% 
  filter(state %in% st,
         result %in% include) %>%
  mutate(popup = paste0("<b>", county, "</b>",
                        "<br/>Population: ", prettyNum(population, big.mark = ","), 
                        "<br/>Hospitals: ", hospitals,
                        "<br/>Expected: ", fit
                        )) %>%  
  mutate(color = case_when(
    result ==  0  ~ at_level,
    result ==  1  ~ over,
    result == -1  ~ under
  )) 

county_map <- counties %>%
  transpose() %>%
  map(~{
    county <- .x
     transpose(county$data)  %>%
       map(~list(
         data = .x$data, 
         color = county$color,
         popup = county$popup
         ))
    }) %>%
  flatten() %>%
  reduce(
    ~ addPolygons(.x, 
                  lng = .y$data$long, 
                  lat = .y$data$lat, 
                  color = .y$color, 
                  popup = .y$popup,
                  weight = 1, fillOpacity = 0.3), 
    .init = initial_map)  


if(add_hospitals == 1) {
  county_map <- county_map %>%
    addCircleMarkers(
      lng = locations$longitude, 
      lat = locations$latitude,
      radius = 1.5 *locations$n,
      popup = locations$popup,
      fillColor = hospital_color, color = "gray", 
      fillOpacity = 0.7,weight = 1)
  } 

county_map

NA

LS0tDQp0aXRsZTogIkFjY2VzcyB0byBjYXJlIg0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KYGBge3J9DQpsaWJyYXJ5KHRpZHl2ZXJzZSkNCmxpYnJhcnkocmxhbmcpDQpsaWJyYXJ5KGxlYWZsZXQpDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbHMgPC0gcmVhZF9jc3YoIi4uL2RhdGEvaG9zcGl0YWxfbGlzdC5jc3YiKQ0KaG9zcGl0YWxzDQpgYGANCg0KYGBge3J9DQpob3NwaXRhbF9sb2NhdGlvbnMgPC0gaG9zcGl0YWxzICAlPiUNCiAgc2VsZWN0KHN0YXRlLCBsb25naXR1ZGUsIGxhdGl0dWRlKQ0Kd3JpdGVfcmRzKGhvc3BpdGFsX2xvY2F0aW9ucywgIi4uL2RhdGEvaG9zcGl0YWxfbG9jYXRpb25zLnJkcyIpDQpob3NwaXRhbF9sb2NhdGlvbnMNCmBgYA0KDQpgYGB7cn0NCnBvcHVsYXRpb24gPC0gdXNtYXA6OmNvdW50eXBvcA0KDQpob3NwaXRhbF9jb3VudCA8LSBob3NwaXRhbHMgJT4lDQogIGNvdW50KHN0YXRlLCBvcmlnaW5hbF9jb3VudHkpIA0KDQoNCnN0YXRlX2hvc3BpdGFscyA8LSBwb3B1bGF0aW9uICU+JQ0KICBsZWZ0X2pvaW4oaG9zcGl0YWxfY291bnQsIGJ5ID0gYygiYWJiciIgPSAic3RhdGUiLCAiY291bnR5IiA9ICJvcmlnaW5hbF9jb3VudHkiKSkgJT4lDQogIG11dGF0ZShob3NwaXRhbHMgPSBpZmVsc2UoaXMubmEobiksIDAsIG4pKSAlPiUNCiAgc2VsZWN0KA0KICAgIGZpcHMsIA0KICAgIHN0YXRlID0gYWJiciwgDQogICAgY291bnR5LCANCiAgICBwb3B1bGF0aW9uID0gcG9wXzIwMTUsDQogICAgaG9zcGl0YWxzDQogICAgKQ0KDQpzdGF0ZV9ob3NwaXRhbHMNCmBgYA0KDQoNCg0KDQpgYGB7cn0NCnNldC5zZWVkKDEwMCkNCg0KaG9zcGl0YWxfc2FtcGxlIDwtIHN0YXRlX2hvc3BpdGFscyAlPiUNCiAgc2FtcGxlX2ZyYWMoMC4zKQ0KDQpob3NwaXRhbF9zYW1wbGUNCmBgYA0KDQpgYGB7cn0NCm1vZGVsIDwtIGxtKGhvc3BpdGFscyB+IHBvcHVsYXRpb24sIGRhdGEgPSBob3NwaXRhbF9zYW1wbGUpDQp3cml0ZV9yZHMobW9kZWwsICIuLi9kYXRhL21vZGVsLnJkcyIpDQpzdW1tYXJ5KG1vZGVsKQ0KYGBgDQoNCmBgYHtyfQ0KcHJlZGljdGlvbnMgPC0gcHJlZGljdChtb2RlbCwgDQogIG5ld2RhdGEgPSBzdGF0ZV9ob3NwaXRhbHMsIA0KICBpbnRlcnZhbCA9ICJwcmVkaWN0aW9uIikgJT4lDQogIGFzX3RpYmJsZSgpICU+JQ0KICBtdXRhdGVfYWxsKHJvdW5kKQ0KDQpwcmVkaWN0aW9ucw0KYGBgDQoNCmBgYHtyfQ0KaG9zcGl0YWxfcmVzdWx0cyA8LSBzdGF0ZV9ob3NwaXRhbHMgJT4lDQogIGJpbmRfY29scyhwcmVkaWN0aW9ucykgJT4lDQogIG11dGF0ZShyZXN1bHQgPSBjYXNlX3doZW4oDQogICAgaG9zcGl0YWxzIDwgbHdyIH4gLTEsDQogICAgaG9zcGl0YWxzID4gdXByIH4gMSwNCiAgICBUUlVFIH4gMCkpICU+JQ0KICBtdXRhdGUoDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eSwgIiBDb3VudHkiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBQYXJpc2giKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiBjaXR5IiksDQogICAgY291bnR5X2tleSA9IHN0cl90b19sb3dlcihjb3VudHlfa2V5KSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlbW92ZV9hbGwoY291bnR5X2tleSwgIiAiKSwNCiAgICBjb3VudHlfa2V5ID0gc3RyX3JlcGxhY2VfYWxsKGNvdW50eV9rZXksICJzdC4gIiwgInNhaW50IiksDQogICAgICApICANCg0Kd3JpdGVfcmRzKGhvc3BpdGFsX3Jlc3VsdHMsICIuLi9kYXRhL2hvc3BpdGFscy5yZHMiKQ0KDQpob3NwaXRhbF9yZXN1bHRzIA0KYGBgDQoNCiMjIFN0YXRlIG1hcA0KDQpgYGB7cn0NCnNoYXBlcyA8LSByZWFkX3JkcygiLi4vZGF0YS9zaGFwZXMucmRzIikgJT4lDQogIG11dGF0ZSgNCiAgICBjb3VudHlfa2V5ID0gc3RyX3RvX2xvd2VyKGNvdW50eSksDQogICAgY291bnR5X2tleSA9IHN0cl9yZW1vdmVfYWxsKGNvdW50eV9rZXksICIgIiksDQogICAgY291bnR5X2tleSA9IHN0cl9yZXBsYWNlX2FsbChjb3VudHlfa2V5LCAic3QuICIsICJzYWludCIpDQogICkgDQoNCnNoYXBlcyANCmBgYA0KYGBge3J9DQpjb3VudHlfaG9zcGl0YWxzIDwtIGhvc3BpdGFsX3Jlc3VsdHMgICU+JQ0KICBzZWxlY3QoLSBjb3VudHkpICU+JQ0KICBsZWZ0X2pvaW4oc2hhcGVzLCBieSA9IGMoInN0YXRlIiwgImNvdW50eV9rZXkiKSkgJT4lDQogIGZpbHRlcihzdGF0ZSAhPSAiUFIiLCBzdGF0ZSAhPSAiQUsiKSAlPiUNCiAgZmlsdGVyKCFpcy5uYShjb3VudHkpKQ0KDQp3cml0ZV9yZHMoY291bnR5X2hvc3BpdGFscywgIi4uL2RhdGEvY291bnR5X2hvc3BpdGFscy5yZHMiLCBjb21wcmVzcyA9ICJneiIpDQoNCmNvdW50eV9ob3NwaXRhbHMNCmBgYA0KDQpgYGB7cn0NCnN0YXRlcyA8LSBjb3VudHlfaG9zcGl0YWxzICU+JQ0KICBncm91cF9ieShzdGF0ZSkgJT4lDQogIHN1bW1hcmlzZSgpICU+JQ0KICBwdWxsKCkNCg0KYGBgDQoNCg0KYGBge3J9DQoNCnVuZGVyIDwtICIjQ0M3OUE3Ig0Kb3ZlciA8LSAiIzAwNzJCMiINCmF0X2xldmVsIDwtICIjMDA4YjAwIg0KaG9zcGl0YWxfY29sb3IgPC0gIiNGMEU0NDIiDQpzdCA8LSBjKCJMQSIpDQojc3QgPC0gc3RhdGVzDQppbmNsdWRlIDwtIGMoLTEsIDEsIDApDQphZGRfaG9zcGl0YWxzIDwtIDANCg0KaW5pdGlhbF9tYXAgPC0gbGVhZmxldCgpICU+JQ0KICBhZGRQcm92aWRlclRpbGVzKHByb3ZpZGVycyRDYXJ0b0RCKSAlPiUNCiAgYWRkTGVnZW5kKCJib3R0b21yaWdodCIsIA0KICAgICAgICAgICAgY29sb3IgID0gYyh1bmRlcixvdmVyICxhdF9sZXZlbCwgaG9zcGl0YWxfY29sb3IpLCANCiAgICAgICAgICAgIGxhYmVscyA9IGMoIkxlc3MgaG9waXRhbHMgdGhhbiBleHBlY3RlZCIsDQogICAgICAgICAgICAgICAgICAgICAgICJNb3JlIGhvc3BpdGFscyB0aGFuIGV4cGVjdGVkIiwgIldpdGhpbiBSYW5nZSIsIA0KICAgICAgICAgICAgICAgICAgICAgICAiSG9zcGl0YWwgTG9jYXRpb24iKSwgDQogICAgICAgICAgICB0aXRsZSAgPSAiTGVnZW5kIixvcGFjaXR5ID0gMC41KQ0KDQpsb2NhdGlvbnMgPC0gaG9zcGl0YWxfbG9jYXRpb25zICU+JQ0KICBmaWx0ZXIoc3RhdGUgJWluJSBzdCkgJT4lDQogIHNlbGVjdChsb25naXR1ZGUsIGxhdGl0dWRlKSAlPiUNCiAgbXV0YXRlX2FsbChyb3VuZCwgMykgJT4lDQogIGNvdW50KGxvbmdpdHVkZSwgbGF0aXR1ZGUpICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIkhvc3BpdGFsczogIiwgbikpDQoNCmNvdW50aWVzIDwtIGNvdW50eV9ob3NwaXRhbHMgJT4lIA0KICBmaWx0ZXIoc3RhdGUgJWluJSBzdCwNCiAgICAgICAgIHJlc3VsdCAlaW4lIGluY2x1ZGUpICU+JQ0KICBtdXRhdGUocG9wdXAgPSBwYXN0ZTAoIjxiPiIsIGNvdW50eSwgIjwvYj4iLA0KICAgICAgICAgICAgICAgICAgICAgICAgIjxici8+UG9wdWxhdGlvbjogIiwgcHJldHR5TnVtKHBvcHVsYXRpb24sIGJpZy5tYXJrID0gIiwiKSwgDQogICAgICAgICAgICAgICAgICAgICAgICAiPGJyLz5Ib3NwaXRhbHM6ICIsIGhvc3BpdGFscywNCiAgICAgICAgICAgICAgICAgICAgICAgICI8YnIvPkV4cGVjdGVkOiAiLCBmaXQNCiAgICAgICAgICAgICAgICAgICAgICAgICkpICU+JSAgDQogIG11dGF0ZShjb2xvciA9IGNhc2Vfd2hlbigNCiAgICByZXN1bHQgPT0gIDAgIH4gYXRfbGV2ZWwsDQogICAgcmVzdWx0ID09ICAxICB+IG92ZXIsDQogICAgcmVzdWx0ID09IC0xICB+IHVuZGVyDQogICkpIA0KDQpjb3VudHlfbWFwIDwtIGNvdW50aWVzICU+JQ0KICB0cmFuc3Bvc2UoKSAlPiUNCiAgbWFwKH57DQogICAgY291bnR5IDwtIC54DQogICAgIHRyYW5zcG9zZShjb3VudHkkZGF0YSkgICU+JQ0KICAgICAgIG1hcCh+bGlzdCgNCiAgICAgICAgIGRhdGEgPSAueCRkYXRhLCANCiAgICAgICAgIGNvbG9yID0gY291bnR5JGNvbG9yLA0KICAgICAgICAgcG9wdXAgPSBjb3VudHkkcG9wdXANCiAgICAgICAgICkpDQogICAgfSkgJT4lDQogIGZsYXR0ZW4oKSAlPiUNCiAgcmVkdWNlKA0KICAgIH4gYWRkUG9seWdvbnMoLngsIA0KICAgICAgICAgICAgICAgICAgbG5nID0gLnkkZGF0YSRsb25nLCANCiAgICAgICAgICAgICAgICAgIGxhdCA9IC55JGRhdGEkbGF0LCANCiAgICAgICAgICAgICAgICAgIGNvbG9yID0gLnkkY29sb3IsIA0KICAgICAgICAgICAgICAgICAgcG9wdXAgPSAueSRwb3B1cCwNCiAgICAgICAgICAgICAgICAgIHdlaWdodCA9IDEsIGZpbGxPcGFjaXR5ID0gMC4zKSwgDQogICAgLmluaXQgPSBpbml0aWFsX21hcCkgIA0KDQoNCmlmKGFkZF9ob3NwaXRhbHMgPT0gMSkgew0KICBjb3VudHlfbWFwIDwtIGNvdW50eV9tYXAgJT4lDQogICAgYWRkQ2lyY2xlTWFya2VycygNCiAgICAgIGxuZyA9IGxvY2F0aW9ucyRsb25naXR1ZGUsIA0KICAgICAgbGF0ID0gbG9jYXRpb25zJGxhdGl0dWRlLA0KICAgICAgcmFkaXVzID0gMS41ICpsb2NhdGlvbnMkbiwNCiAgICAgIHBvcHVwID0gbG9jYXRpb25zJHBvcHVwLA0KICAgICAgZmlsbENvbG9yID0gaG9zcGl0YWxfY29sb3IsIGNvbG9yID0gImdyYXkiLCANCiAgICAgIGZpbGxPcGFjaXR5ID0gMC43LHdlaWdodCA9IDEpDQogIH0gDQoNCmNvdW50eV9tYXANCiAgDQpgYGANCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg0K